


#################################################
#### Analysis for State Senate Districts ########
#################################################


# District id variable
districts<-read.csv(paste("districts_ssd.csv",sep=""))
district.uid<-districts$ssd_district_uid


### Prepare survey data for response model with state senate variables 
data1<-read.dta("gay_marriage_senate.dta")
state.gaymarriage<-data1$state_gaymarriage
race2<-data1$race
cst2<-factor(data1$cst)
education5<-factor(data1$education)
factor.education<-factor(data1$education5)
region<-data1$region_num
districturban.full<-data1$ssd_district_urban
districtincome.full<-data1$ssd_district_income/100000
districtveteran.full<-data1$ssd_district_veteran
statereligion.full<-data1$state_religion
stateunion.full<-data1$state_union
districtgay.full<-data1$ssd_district_gay
district_num<-data1$ssd_district_uid
gender<-data1$gender


age <- data1$age_exact






age.group <- rep(1,length(age))
age.group[age<20] <- 1
age.group[age>19&age<25] <- 2
age.group[age>24&age<30] <- 3
age.group[age>29&age<35] <- 4
age.group[age>34&age<40] <- 5
age.group[age>39&age<45] <- 6
age.group[age>44&age<50] <- 7
age.group[age>49&age<55] <- 8
age.group[age>54&age<60] <- 9
age.group[age>59&age<65] <- 10
age.group[age>64&age<70] <- 11
age.group[age>69&age<75] <- 12
age.group[age>74&age<80] <- 13
age.group[age>79&age<85] <- 14
age.group[age>84&age<90] <- 15
age.group[age>89&age<95] <- 16
age.group[age>94] <- 17




### Prepare data for predictions with state senate data (n=77680, which is 2*4*5 * 1942) 
PS.Data<-read.csv(paste("Post-Stratification Data - RImport (ssd).csv",sep=""))
proportion<-PS.Data$proportion
ps_cregion<-PS.Data$ps_cregion
cst_ps<-PS.Data$cst_ps
cat_educ<-PS.Data$cat_educ
cat_race<-PS.Data$cat_race
cat_female<-PS.Data$cat_female
ps_district_income <-PS.Data$ps_district_income/100
ps_district_urban<-PS.Data$ps_district_urban
ps_district_veteran<-PS.Data$ps_district_veteran
ps_districtgay<-PS.Data$ps_district_gay
ps_state_religion<-PS.Data$ps_state_religion
ps_state_union<-PS.Data$ps_state_union
ps_district_uid<-PS.Data$ssd_district_uid


# Prepare data for disaggregation (calculate mean for each senate district)
sample.mean.gaymarriage<-tapply(state.gaymarriage, data1$ssd_district_uid, mean)
sample.count.gaymarriage<-tapply(state.gaymarriage, data1$ssd_district_uid, length)
sample.mean.gaymarriage.names<-as.numeric(rownames(sample.mean.gaymarriage))
sample.mean.gaymarriage.adj1<-cbind(sample.mean.gaymarriage, sample.mean.gaymarriage.names)
sample.mean.gaymarriage.adj2<-data.frame(matrix(NA, nrow=1942, ncol=4))

sample.count.gaymarriage.names<-as.numeric(rownames(sample.count.gaymarriage))
sample.count.gaymarriage.adj1<-cbind(sample.count.gaymarriage, sample.count.gaymarriage.names)
# 1 senate district is missing


## Create a vector of district effects without any missing cells (zeros where no data for a district)
for (i in 1:length(unique(district.uid)))
{
sample.mean.gaymarriage.adj2[i,1]<-sum(sample.mean.gaymarriage.adj1[sample.mean.gaymarriage.names==i,1], na.rm=FALSE)
sample.mean.gaymarriage.adj2[i,2]<-sum(sample.count.gaymarriage.adj1[sample.count.gaymarriage.names==i,1], na.rm=FALSE)
}

sample.mean.gaymarriage.adj2[,3]<-districts$cst

sample.mean.gaymarriage.adj2[,4]<-sample.mean.gaymarriage.adj2[,1]

for (i in 1:length(unique(district.uid)))
{
if (sample.mean.gaymarriage.adj2[i,2]==0){
sample.mean.gaymarriage.adj2[i,4]<-NA
}
}
sample.mean.gaymarriage<-sample.mean.gaymarriage.adj2[,4]


## Marriage response model

DATAmodel <- data.frame(cbind(age.group, age, state.gaymarriage, race2, gender, cst2, district_num, education5, districtincome.full, districturban.full, districtveteran.full, statereligion.full, stateunion.full, region, districtgay.full))

#marriage.model1 <-  glmer(formula = state.gaymarriage ~ (1 | race2) + gender + (1 | cst2) + (1 | district_num) +  (1 | education5)+ districtincome.full + districturban.full + districtveteran.full + statereligion.full + stateunion.full + (1 | region) + districtgay.full, family=binomial(link="logit"))
marriage.model <-  glmer(formula = state.gaymarriage ~ (1 | race2) + gender + (1 | cst2) + (1 | district_num) +  (1 | education5)+ districtincome.full + districturban.full + districtveteran.full + statereligion.full + stateunion.full + (1 | region) + districtgay.full, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa"))

# Get RE for each district (1 is missing)
district.effects<-ranef(marriage.model)$district_num
district.effects.names<-as.numeric(rownames(district.effects))
district.effects2<-cbind(district.effects, district.effects.names)

district.effects3<-data.frame(matrix(NA, nrow=1942, ncol=1))

# Create a vectors of district effects without any missing cells (zeros where no data for a district)
for (i in 1:length(unique(district.uid))){
	j<-"i"
	district.effects3[i,1]<-sum(district.effects2[district.effects2$district.effects.names==i,1], na.rm=TRUE)
}


# Get RE for each state 
state.effects<-ranef(marriage.model)$cst2
state.effects.names<-rownames(state.effects)
state.effects2<-cbind(state.effects, state.effects.names)
state.effects3<-data.frame(matrix(NA, nrow=51, ncol=2))
states<-unique(cst_ps)
state.effects3[,2]<-states

state.effects3[1:48,1]<-state.effects2[,1]
states2<-factor(states)

## Create a vectors of district effects without any missing cells (zeros where no data for a district)
for (i in 1:length(states)){
	state.effects3[i,1]<-sum(state.effects2[state.effects.names==states[i],1], na.rm=TRUE)
}
state.effects3[50,2]<-"AK"
state.effects3[51,2]<-"HI"

# Make predictions
cell.pred<-invlogit(fixef(marriage.model)["(Intercept)"]+
ranef(marriage.model)$education5[cat_educ,1] + ranef(marriage.model)$race2[cat_race,1] +   fixef(marriage.model)["gender"] *cat_female +    district.effects3[ps_district_uid,1] + state.effects3[cst_ps,1]  +  ranef(marriage.model)$region[ps_cregion,1]  + fixef(marriage.model)["districtincome.full"]*ps_district_income + fixef(marriage.model)["districturban.full"]*ps_district_urban + fixef(marriage.model)["districtveteran.full"]*ps_district_veteran+fixef(marriage.model)["statereligion.full"]*ps_state_religion + fixef(marriage.model)["stateunion.full"]*ps_state_union+ fixef(marriage.model)["districtgay.full"]*ps_districtgay)


# Post-stratification
cell.pred.weighted <- cell.pred * proportion

mrp.mean.gaymarriage<-tapply(cell.pred.weighted, PS.Data$ssd_district_uid, sum)
write.csv(mrp.mean.gaymarriage,file="stategaymarriage.csv")

mrp.mean.gaymarriage2<-data.frame(matrix(NA, nrow=1942, ncol=3))
mrp.mean.gaymarriage2<-cbind(mrp.mean.gaymarriage, districts$ssd_fips_num)



# Get disaggregation and MrP data for the congressional districts for each of the 5 states as well as the referendum results from the datasets

ohio <-read.csv(paste("OH_GM_SSD.csv",sep=""))
mrp.ohio<-mrp.mean.gaymarriage2[1389:1421,]
raw.ohio<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.ohio<-sample.mean.gaymarriage[1389:1421]

ca <-read.csv(paste("CA_GM_SSD.csv",sep=""))
mrp.ca<-mrp.mean.gaymarriage2[121:160,]
raw.ca<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.ca<-sample.mean.gaymarriage[121:160]

wi<-read.csv(paste("WI_GM_SSD.csv",sep=""))
mrp.wi<-mrp.mean.gaymarriage2[1880:1912,]
raw.wi<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.wi<-sample.mean.gaymarriage[1880:1912]

mi <-read.csv(paste("MI_GM_SSD.csv",sep=""))
mrp.mi<-mrp.mean.gaymarriage2[815:852,]
raw.mi<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.mi<-sample.mean.gaymarriage[815:852]

az <-read.csv(paste("AZ_GM_SSD.csv",sep=""))
mrp.az<-mrp.mean.gaymarriage2[56:85,]
raw.az<-data.frame(matrix(NA, nrow=33, ncol=1))
raw.az<-sample.mean.gaymarriage[56:85]
#raw.az[15]<-NA


ssd.raw<-append(raw.ohio, raw.ca)
ssd.raw<-append(ssd.raw, raw.wi)
ssd.raw<-append(ssd.raw, raw.mi)
ssd.raw<-append(ssd.raw, raw.az)
ssd.mrp<-append(mrp.ohio[,1], mrp.ca[,1])
ssd.mrp<-append(ssd.mrp, mrp.wi[,1])
ssd.mrp<-append(ssd.mrp, mrp.mi[,1])
ssd.mrp<-append(ssd.mrp, mrp.az[,1])
ssd.referendum<-append(ohio[,2], ca[,2])
ssd.referendum<-append(ssd.referendum, wi[,2])
ssd.referendum<-append(ssd.referendum, mi[,2])
ssd.referendum<-append(ssd.referendum, az[,2])
ssd.regression1<-lm(ssd.referendum~ssd.raw)
ssd.regression2<-lm(ssd.referendum~ssd.mrp)


se.dis <- (ssd.raw - ssd.referendum)^2
se.mrp <- (ssd.mrp - ssd.referendum)^2


